Synchronised firing induced by network dynamics in excitable 
systems 
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Abstract -We study the collective dynamics of an ensemble of coupled identical FitzHugh-Nagumo elements 
in their excitable regime. We show that collective firing, where all the elements perform their individual 
firing cycle synchronously, can be induced by random changes in the interaction pattern. Specifically, on a 
sparse evolving network where, at any time, each element is connected with at most one partner, collective 
firing occurs for intermediate values of the rewiring frequency. Thus, network dynamics can replace noise 
and connectivity in inducing this kind of self-organised behaviour in highly disconnected systems which, 
otherwise, wouldn't allow for the spreading of coherent evolution. 



Introduction. - The spontaneous emergence of coherent 
behaviour in populations of interacting units — be they of phys- 
ical, chemical, biological, or technological nature — is crucial 
to their collective function. Synchronisation of several kinds 
occurs in such disparate systems as mechanical oscillators, 
bio-molecular reactions, neural networks, insect societies, hor- 
monal cycles, coupled lasers, and Josephson junctions. Math- 
ematical models for this variegated class of phenomena have 
been proposed in terms of ensembles of coupled dynamical sys- 
tems of different types: linear and nonlinear oscillators, chaotic 
elements, excitable units, among others 0121. 

It is a well-established fact that, in a population of interact- 
ing elements, sufficiently strong, attractive coupling induces 
self-organised synchronisation. This occurs even in the pres- 
ence of external noise, or when the individual behaviour of 
each element is chaotic, or when elements are not identical 
to each other, with the proviso that the population is well- 
interconnected in such a way that information about the state 
of any element can reach any other. While the structure of the 
interaction pattern can affect details in the collective dynamics 
|3,4|, connectedness and strong coupling generally guarantee 
synchronisation. 

It has recently been shown that, both in synchronisation and 
in contact processes (such as epidemics spreading), instanta- 
neous lack of connectivity can be compensated by dynamical 
rewiring of the interaction network J5]|6l. Specifically, in pop- 
ulations with very sparse, disconnected instantaneous interac- 
tion patterns, the respective transitions to full synchronisation 



and to endemic states are triggered by increasing the rewiring 
rate. This result is relevant, especially, to biological and social 
networks, where potential contacts between the members of a 
population are not continuously realised, but can occasionally 
be activated. 

In this Letter, we disclose a related but different phe- 
nomenon, concerning the collective dynamics of populations 
of excitable units on evolving networks. Interacting excitable 
elements, which individually perform a "firing" cycle in phase 
space if perturbed strongly enough from their quiescent state, 
are known to undergo collective synchronised firing induced 
by external noise QUO and by repulsive interactions (3J- We 
show here that, even in the absence of noise, an ensemble of 
coupled FitzHugh-Nagumo excitable elements on an evolving, 
very sparse, network exhibits collective firing for intermedi- 
ate values of the rewiring rate. This phenomenon in charac- 
terised numerically for different network topologies and semi- 
quantitatively explained in terms of the perturbations that net- 
work reconnections impose on the individual dynamics of each 
element. 

Excitable elements on an evolving network. - FitzHugh- 
Nagumo excitable elements constitute an archetypical model 
for type-II excitability, which occurs in many natural and arti- 
ficial systems ranging from epidemic spreading, to neural and 
cardiac tissues 110111 11 to chemical reactions and electronic de- 
vices 1121 . The model is defined in terms of an activatory (fast) 
variable x and an inhibitory (slow) variable y. We consider 
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an ensemble of yV sparsely connected FitzHugh-Nagumo ele- 
ments, whose dynamics is given by 



isfies 



yj 



(1) 

(2) 



for j = 1 , . . . , N. The time-dependent factor Kj(t) weights the 
interaction between elements j and f, as explained below. The 
small parameter e measures the time-scale ratio between the 
fast and slow variables Xj and yj. The positive parameter a 
characterises the dynamical regime of the non-interacting (K = 
0) element: for a < 1, it performs a periodic oscillation in the 
(xj,yj)-plane, while for a > 1 its behaviour is excitable. In 
this latter regime, and in the absence of external perturbations, 
the non-interacting element asymptotically approaches the sole 
stable fixed point (x eq ,y eq ) = (a,a 3 /3-a) and remains quiescent 
there. Under a sizeable perturbation, however, the element may 
exit the vicinity of the fixed point and return to it after a long 
excursion in phase space — usually referred to as a firing cycle, 
or spike. 

Our FitzHugh-Nagumo elements interact through a sparse 
evolving network such that, at any given time, each element 
is coupled to at most one partner. The coupling constant in 
eq. (Q]i is Kj(t) = k when element j interacts with a generic 
partner j* (not necessarily the same at all times), and Kj(t) = 
when j is isolated. This rather extreme sparseness determines, 
in a sense, the most unfavourable situation for the emergence of 
collective phenomena in an ensemble of interacting units. We 
expect network dynamics to replace connectivity in triggering 
collective behaviour. 

In our numerical simulations, we study two different 
schemes for the network dynamics. In the first one, the net- 
work consists of exactly N/2 undirected links, distributed in 
such a way that every element is always coupled to exactly one 
partner. As time elapses, two connected pairs of elements are 
occasionally chosen at random to mutually exchange their part- 
ners. Thus, two links in the network are rewired. 

The second scheme for network dynamics is built on top of 
an underlying (undirected, connected) network G with a fixed 
number of links. During the dynamical evolution, however, 
only a subset of the links is active. The links of the underly- 
ing network G thus represent the potential connections in the 
actual interaction pattern. The initial network is generated by 
successively selecting elements in a random order. If the cho- 
sen element is isolated, the link to one of its still isolated neigh- 
bours in G gets activated. If no available neighbours exist, the 
element remains with no active connection. During evolution, 
an inactive link from G is occasionally chosen at random and 
gets activated. At the same time, pre-existing links of the newly 
connected elements become deactivated. 

It is not difficult to realise that, if there is no correlation be- 
tween the degrees of neighbour nodes in the underlying net- 
work G, the frequency a>+(z) with which an isolated node of de- 
gree z becomes connected, exactly equals the frequency o»_(z) 
with which it becomes isolated when it is connected. In turn, 
the probability P to find the node connected to any partner sat- 



P = w+(z)(l - P) - co-(z)P- 



(3) 



Therefore, for asymptotically long times, P — 1/2 for any z. In 
other words, in our second reconnection scheme and for long 
times, there are on the average N/2 connected elements — and, 
consequently, N/4- links — at any time. The resulting interac- 
tion pattern is thus twice as sparse as in the first scheme. 

In both reconnection schemes, we denote by A the reconnec- 
tion rate, i.e. the probability per time unit that the partner of 
any given element changes. 

Order parameters. - To characterise the collective prop- 
erties of our system in the framework of the standard theory 
of synchronised oscillators |[T3l . it is convenient to compute a 
quantity describing the phase of each excitable element along 
its firing cycle. For the model defined by eqs. (JTJ and (fJJ, 
the excursion in phase space occurs around the origin of the 
(x 7 -,y ; )-plane. Thus, a suitable definition of the phase is simply 



<pj(t) = tan 



yj(t) 



Xj(t) 



(4) 



The behaviour of the ensemble, including possible transi- 
tions between different collective dynamical regimes, can be 
statistically characterised by a pair of order parameters defined 
in terms of the individual phases <pj(t) [8 9 |. First, we take the 
average of the location of the particles on the unit circle, 
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p(f)exp[/ v F(f)] = - 2j exp[/0;(f)] 



(5) 
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and compute the Kuramoto order parameter as p = (p(t)), 
where (•) stands for the time average over a long time inter- 
val [TT~3| . This parameter measures the degree of synchroni- 
sation attained by the ensemble: with full synchronisation we 
have p = 1, whereas for a state where phases are uniformly 
distributed over [0, 2k) we have p ~ N~ 1 ^ 2 . 

In excitable systems, however, the Kuramoto order parame- 
ter does not allow to discern between the case where phases are 
statically synchronised at the fixed point <p eq = tan~ l (y eq /x eq ), 
and the case where they rotate coherently, as expected to occur 
in the regime of collective firing. To discriminate between static 
and dynamic synchronisation, we compute the Shinomoto- 
Kuramoto order parameter |[T4l . 

£ = (|p(/) exp[/T(f)] - (p(t) «p[iT(0]>|) , (6) 

which differs from zero for synchronous firing only. 

A third relevant order parameter, frequently used in the anal- 
ysis of stochastic transport |fT51 , is the current, which we com- 
pute as 



./ 



(if 



Xj(t) 



(7) 



i.e. as the time average of the absolute mean velocity along the 
coordinate x. It gives a measure of the level of (not necessarily 
synchronised) firing in the ensemble. 
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Fig. 1: (Colour on-line) The order parameters as functions of the 
reconnection rate A, for the first reconnection scheme, with k = 1, 
a = 1.02, and e = 1CT 3 . Upper panel: the Kuramoto order parameter 
p; central panel: the Shinomoto-Kuramoto order parameter f ; Lower 
panel: the current J. Different curves correspond to different system 
sizes: N = 10 (O), 50 (□), 200 (0), 800 (A), and 2000 (<). Joining 
lines are plotted as a guide to the eye. 



Numerical results. - We have performed extensive com- 
puter simulations of the model defined by eqs. (Q} and (O for 
a - 1.02 (excitable regime) and k = 1, with the corresponding 
network dynamics. The results presented here are qualitatively 
representative of a broad parameter range in the same regime. 
Order parameters were computed after the system reached a 
stationary state, with no further changes in its dynamical be- 
haviour. 

In fig. [1] we show numerical results for the order parameters 
as functions of the reconnection rate A, for the first reconnection 
scheme and various system sizes. Both for small and large val- 
ues of A, the Kuramoto order parameter is equal to one, while 
the Shinomoto-Kuramoto order parameter and the current van- 
ish. This situation corresponds to a state where all the elements 
of the ensemble are at rest at the same point in phase space, 
namely, the fixed point (x eq ,;y e q). 

For intermediate reconnection rates, on the other hand, we 
find an interval where the Kuramoto order parameter p < 1, 
indicating that the elements are distributed over phase space. 
Within the same interval, both the Shinomoto-Kuramoto or- 
der parameter £ and the current J become positive and attain 
considerably high maxima. This is an indication of collective 
firing with a concurrent phase-space flow, and constitutes our 
main finding: reconnection events at intermediate rates induce 
self-organised coherent behaviour in an otherwise disconnected 
ensemble of FitzHugh-Nagumo excitable elements. 

Figure Q] also shows that the order parameters become inde- 
pendent of the system size as N grows. This suggests that the 
regime of collective firing exists even in the thermodynamic 
limit. 

Figure|2]displays the order parameters for different values of 
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Fig. 2: (Colour on-line) As in fig.[T| f° r a system of size N = 400 and 
different values of the time-scale ratio: e = 1 (O), 10 _1 (□), 10~ 2 (0), 
10~ 3 (A), 10~ 4 (<), and 10~ 5 (V). 



the time-scale ratio e in a system of size N = 400. These re- 
sults show that, when there is no difference in the time scales 
associated to the variables x and y (e = 1), collective firing is 
absent and all the elements remain quiescent at the fixed point. 
As e decreases, however, the phenomenon takes place for inter- 
mediate values of A and seems to approach a well-defined limit 
for € — > 0. We provide a semi-quantitative analysis of this limit 
in the next section. 

To relax the condition that every element has a partner at 
any time, we have used our second reconnection scheme with 
two kinds of topologies for the underlying network G. We re- 
call that, with this scheme, the resulting instantaneous interac- 
tion pattern is more sparse than in the previous case. Firstly, 
we have considered a scale-free underlying network generated 
by the preferential attachment rule |fl6l , where each added 
node is connected to m pre-existing nodes. Secondly, we have 
taken a small-world network built up from the rewiring, with 
probability p, of the links of a two-dimensional network with 
Moore neighbourhood, following the Watts-Strogatz prescrip- 
tion ifTTll . We have numerically verified that, as advanced 
above, the average number of connected elements at long times 
fluctuates around N/2. For the scale-free networks, this number 
is slightly, but systematically, larger, which can be attributed 
to spurious degree-degree correlations in the highly heteroge- 
neous degree distribution generated by preferential attachment 
in our finite-size system. 

Figure [3] shows the order parameters for a system of size 
N = 400, underlain by scale-free networks with two values 
of m and small-world networks with two values of p. Solid 
curves stand for the corresponding results for the first recon- 
nection scheme. Overall, the results are largely independent of 
both the reconnection scheme and the topology of the under- 
lying network and, consequently, of the number of connected 
elements. 
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Fig. 3: (Colour on-line) As in fig. Q] for the second reconnection 
scheme, with e = 10~ 3 and N = 400. Open symbols: scale-free un- 
derlying network generated by preferential attachment with m = 3 (<>) 
and 10 (□). Full symbols: small-world underlying networks generated 
by rewiring with p = 0.01 (•) and 0.1 (■). The solid line corresponds 
to the results for the first reconnection scheme. 



Interpretation. - Whereas a full analytical description 
of synchronised firing in dynamical networks of FitzHugh- 
Nagumo excitable elements seems to be out of reach, it is possi- 
ble to sketch a semi-quantitative picture that plausibly explains 
the occurrence of this collective phenomenon for intermedi- 
ate values of the network reconnection rate. The following 
arguments focus on our first reconnection scheme, but can be 
straightforwardly extended to the second. 

Consider eqs. (fTJ and (|2]) for e — > 0. In this limit, and in the 
absence of interactions (k = 0), the fast variable Xj(f) follows 
adiabatically the slow variable yj(t) along the nullcline xj = 0. 
Let us introduce the auxiliary variable 



1 , 

T]j(Xj) = -A , - Xj 



(8) 



which, along the nullcline, satisfies 77 y = yj + k\Xf - Xjj. For 
k = 0, we have yj = r\j. Differentiating with respect to time and 
taking into account eq. (O yields 



(9) 



with £ = £f - xj and Xj(T]j) given by the inverse of the function 
in eq. ([SJ. Note that Xj{rjj) is defined piecewise, depending on 
how Xj compares with ±1. 

The arrowed bold lines in fig.[4]represent the phase-space tra- 
jectories of an non-interacting element. In the limit e — » 0, it is 
always found on the stable branches (either Xj < -1 or xj > 1) 
of the nullcline, and asymptotically approaches the fixed point 
at (x e q,??eq) = (a,a 3 /3 - a), plotted in the figure as an empty 
dot. If, as illustrated by the grey arrow, the element is perturbed 
from the fixed point toward negative values of r\j and beyond 
the minimum r\j{xj — 1) = -2/3, it immediately reaches the 
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Fig. 4: Phase-space dynamics of a single non-interacting FitzHugh- 
Nagumo excitable element j on the (xj,Tjj) plane. The narrow line 
represents eq. Arrowed bold lines are possible trajectories of the 
element. The fixed point (x eq , 77 eq ) is represented as an empty dot. The 
grey arrow schematises the effect of a perturbation from the fixed point 
toward negative values of rjj. 



leftmost stable branch and begins its excursion upwards. When 
it reaches rjj(xj = -1) = 2/3, it jumps to the rightmost branch 
and, from then on, it moves toward the fixed point. The firing 
cycle is thus completed. Integration of eq. (0 with k = and 
a > 1 shows that, if the leftmost branch is reached at rjj as -2/3, 
the time spent on that branch is Ti e ft = 1/2 + 0[a - 1]. In turn, 
the typical time for relaxation toward the fixed point on the 
rightmost branch is Tri e ht - a 2 - I. 

Consider now the effect of interaction (k + 0) on the in- 
dividual dynamics of element j. If the reconnection rate A is 
sufficiently small, so that /T 1 » Ti e ft, Tnght! element j remains 
connected to the same partner j* over times which are long as 
compared with the typical time scales needed to reach the vicin- 
ity of (x eq , ?7eq). Irrespectively of the value of k, the two coupled 
elements approach the fixed point well before their mutual link 
breaks and they are reconnected to different partners. When re- 
connection finally happens, however, all elements will be found 
near the fixed point and the change of partner will have essen- 
tially no effect on the subsequent dynamics of j. Therefore, 
the whole ensemble converges to (jc e q, i]eq) over times of order 
"Heft + T rfght and remains there indefinitely. For sufficiently small 
A, hence, sustained collective firing is absent. 

As A grows and reconnection becomes more frequent, the 
term k^j in the right-hand side of eq. (0 acquires the charac- 
ter of a fluctuating force, analogous to additive noise. Since 
^j(t) = Xj*(t) - Xj(t), its time dependence consist of a relatively 
smooth variation along the periods where element j's partner j* 
does not change, punctuated by sharp delta-like "kicks" when 
reconnection occurs. Even if j has already reached the vicinity 
of (jt e q, rjeq), a kick due to reconnection with an element which 
is transiting the leftmost branch region may force /' to move 
away from the fixed point and reinitiate its firing cycle. This 
event is schematised by the grey arrow in fig. [4] At appropri- 
ate values of the reconnection rate, with most of the ensem- 
ble near (x eq , j] eq ), just a few "outliers" along the firing cycle 
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Fig. 5: (Colour on-line) Upper panel: Standard deviation cr f of the 
"noise" gj(f) as a function of the reconnection rate A, for k = 1, a = 
1.02, A' = 400, and different values of the of the time-scale ratio: 
e = 1 (O), 10-' (□), 10~ 2 (0), 10~ 3 (A), 10~ 4 (<), and 10~ 5 (V). Lower 
panels: Time evolution of the coordinate xj and the "noise" of three 
randomly selected elements, for e = 10~ 3 , with A = 0.32 (upper row) 
and A = 3.2 (lower row). 



are able to induce a cascade of transitions from the fixed point 
to the cycle, and collective firing is thus triggered. Our nu- 
merical results show that, precisely, collective firing occurs for 

A> \ ~ (Tieft + Tright)" 1 . 

If reconnection grows even more frequent, within the time 
scales relevant to the dynamics of a single element, the "noise" 
term kgj averages out to its mean value over the whose ensem- 
ble. Therefore, each element is effectively subject to the action 
of the average state of the ensemble. In this situation, the in- 
teraction between elements is equivalent to global (all-to-all) 
coupling [5 ,6|. Since all the elements are identical, global cou- 
pling leads the ensemble to collapse to the fixed point J2J, and 
collective firing is thus suppressed. 

The upper panel of fig. [5] shows the standard deviation cr^ 
of the "noise" £-j(t), averaged over the ensemble and over time, 
as a function of the reconnection rate A and for various val- 
ues of the time-scale ratio e. The lower panels show the co- 
ordinate xj(t) and the "noise" £ 7 (f) as functions of time for a 
few selected elements, and two values of the reconnection rate: 
A = 0.32 (upper row), which corresponds to the threshold of 
global firing, and X = 3.2 (lower row), where global firing is 
well developed. For the latter, the synchronous pulsing of the 
coordinate xj(t) is apparent. 

Conclusion. - Synchronised collective firing in ensembles 
of coupled excitable elements was known to be triggered by 
external noise QUO and by disorder in the interaction pattern 
||9l — in this latter case, due to the simultaneous presence of at- 
tractive and repulsive interactions. In both situations, the emer- 
gence of this form of collective behaviour requires the intensity 
of noise or the degree of disorder to be neither too small nor too 



high: it is at an intermediate level of fluctuations that the sys- 
tem has the appropriate dynamical flexibility as to self-organise 
into coherent evolution. 

In this Letter, we have shown that, in the absence of ex- 
ternal noise, the fluctuations associated with network dynam- 
ics — when the interaction pattern is rewired with a certain 
frequency — are as well able to induce collective firing of cou- 
pled excitable elements. As in the previous instances, coherent 
evolution is observed for intermediate values of the rewiring 
frequency. In the present situation, network dynamics has the 
crucial additional role of replacing the connectivity necessary 
to warrant the spreading of information about the individual 
states of the excitable elements all over the ensemble. In fact, 
by construction, the instantaneous interaction pattern is highly 
diluted, with one or less neighbour connected to each element 
at any time. This effect of network dynamics had already been 
pointed out in chaotic synchronisation and in contact processes 
[5 6|. Our results suggest that the phenomenon of collective 
firing is remarkably independent of the underlying structure of 
interactions and has a well-defined behaviour in the thermody- 
namic limit of infinitely large ensembles. 

The present analysis pertains to the study of the ample va- 
riety of systems where interaction patterns are not static, but 
change with time either driven by external influences or in re- 
sponse to the state of the system itself, or as a combination of 
both effects. While this important dynamical aspect of com- 
plex systems has often been disregarded, our results — among 
other recent work — highlight its role in the emergence of self- 
organised collective evolution. 
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